Cube with vibrating sides (Interior)

using LinearAlgebra, BoundaryIntegralEquations, IterativeSolvers, Plots, SpecialFunctions, Meshes
import WGLMakie as wgl # WGLMakie integrates into VSCode. Other backends can also be used.
wgl.set_theme!(resolution=(800, 800))

Loading and visualizing the triangular cube mesh

First we define the path to the mesh file

mesh_file = joinpath(dirname(pathof(BoundaryIntegralEquations)),"..","examples","meshes","1m_cube_extremely_coarse");

Now we read in the mesh

physics_orders  = [:linear,:geometry,:disctriconstant,:disctrilinear,:disctriquadratic];
mesh = load3dTriangularComsolMesh(mesh_file;geometry_order=:linear, physics_order=physics_orders[5])

Number of elements: 	 84
Number of unkowns:  	 504
Geometry defined by:	 TriangularLinear{Float64}
Physics defined by: 	 DiscontinuousTriangularQuadratic{Float64}

We now define the mesh entity that contain the boundary condition. In this case it is boundary 0.

bc_ents = [0];

Now two simple meshes is created. One for the boundary condition and one for the remaining part of the mesh

simple_bc   = create_bc_simple_mesh(mesh,bc_ents);
simple_mesh = create_bc_simple_mesh(mesh,bc_ents,false);

Using the simple meshes we can visualize the mesh, with boundary 0 (where velocity condition will be applied) shown in red

viz(simple_mesh;showfacets=true)
viz!(simple_bc;showfacets=true,color=:red)
wgl.current_figure()

Setting up constants

frequency = 100.0;      # Frequency                [Hz]
c  = 343.0;             # Speed up sound           [m/s]
ρ₀ = 1.21;              # Mean density             [kg/m^3]
Z₀ = ρ₀*c;              # Characteristic impedance [Rayl]
v₀ = 1.0;               # Speed in the x-direction [m/s]
k  = 2π*frequency/c;    # Wavenumber

Analytical Solution

The analytical description of the interior pressure in unit cube with the side at $x=0$ be applied a velocity of $v_{0}$

\[ p_\text{analytical}(x) = -\frac{\mathrm{i}Z_0v_{0}\cos(k(1 - x))}{\sin(k)}\]

where $Z_0$ is the characteristic impedance and $k$ is the wavenumber. We now compute the analytical expression is computed at the points

x_ana = collect(0.00:0.01:1)
p_analytical = -im*Z₀*v₀*cos.(k*(1 .- x_ana))/(sin(k));

Solution using the dense BEM

We start by solving the BEM system using dense matrices. For this we need to first assemble the matrices

F,G,C = assemble_parallel!(mesh,k,mesh.sources,n=2,m=2,progress=false);
H  = Diagonal(C) - F; # Adding the integral free term

In order to apply the boundary coniditon we must find the nodes corresponding to surface entity 0

bc_elements = 0 .== mesh.entities;
bc2 = sort(unique(mesh.physics_topology[:,bc_elements]));

We now define a vector corresponding to the velocities at bc2

v_bem = zeros(ComplexF64,length(C));
v_bem[bc2] .= im*Z₀*k*v₀*mesh.normals[1,bc2];

Using this we can define the right-hand side

b = G*v_bem;

Finally the pressures can be computed by solving a linear system of equations

p_bem = gmres(H,b;verbose=true);

=== gmres ===
rest	iter	resnorm
  1	  1	1.87e+03
  1	  2	9.34e+02
  1	  3	5.08e+02
  1	  4	1.82e+02
  1	  5	8.10e+01
  1	  6	3.87e+01
  1	  7	1.92e+01
  1	  8	7.68e+00
  1	  9	3.56e+00
  1	 10	1.37e+00
  1	 11	6.31e-01
  1	 12	2.55e-01
  1	 13	1.52e-01
  1	 14	6.56e-02
  1	 15	4.23e-02
  1	 16	1.67e-02
  1	 17	8.77e-03
  1	 18	3.87e-03
  1	 19	1.79e-03
  1	 20	1.09e-03
  2	  1	7.55e-04
  2	  2	4.27e-04
  2	  3	3.52e-04
  2	  4	3.14e-04
  2	  5	1.64e-04
  2	  6	8.69e-05
  2	  7	3.65e-05
  2	  8	1.99e-05

Solution using the FMM-BEM

Similarly we can compute the BEM solution using the Fast Multipole Operators

Gf = FMMGOperator(mesh,k,depth=2);
Ff = FMMFOperator(mesh,k,depth=2);
Hf = 0.5I - Ff;
bf = Gf*v_bem;                     # Computing right-hand side
p_fmm = gmres(Hf,bf;verbose=true); # Solving the linear system

=== gmres ===
rest	iter	resnorm
  1	  1	1.87e+03
  1	  2	7.28e+02
  1	  3	1.83e+02
  1	  4	5.64e+01
  1	  5	1.38e+01
  1	  6	2.95e+00
  1	  7	9.62e-01
  1	  8	2.50e-01
  1	  9	6.49e-02
  1	 10	2.44e-02
  1	 11	7.74e-03
  1	 12	2.61e-03
  1	 13	1.01e-03
  1	 14	2.21e-04
  1	 15	7.88e-05
  1	 16	2.60e-05

Field evaluations

In order to compare the results with the analytical solution we must use the found surface pressures to compute pressure at the interior field points. We therefore start by creating a matrix of points (as columns). Below $N$ linearly spaced points defined by $(x,0.5,0.5)$ with $x\in[0.1,0.9]$ is created.

N = 20;
x = collect(0.1:(0.9-0.1)/(N-1):0.9);
y = 0.5ones(N);
z = 0.5ones(N);
X_fieldpoints = [x'; y'; z'];

Using the described field-points we can assemble the (dense) matrices for the field point evaluation

Fp,Gp,Cp = assemble_parallel!(mesh,k,X_fieldpoints,n=5,m=5,progress=false);

Using the matrices we can evalute the pressure at the field points using the previously found surface pressure

p_field  = Fp*p_bem + Gp*v_bem;
pf_field = Fp*p_fmm + Gp*v_bem;

Plotting the field point pressure it can be seen that all 3 methods find the correct pressures

plot(x_ana,abs.(p_analytical./Z₀), label="Analytical")
scatter!(x,abs.(p_field./Z₀),  label="Full", markershape=:rect)
scatter!(x,abs.(pf_field./Z₀), label="FMM")
ylabel!("|p/Z₀|"); xlabel!("r/a")

The surface pressures can also be plotted. Note that it is only possible to plot linear meshes, meaing that we must remove the quadratic parts.

data_mesh,data_viz = create_vizualization_data(mesh,p_fmm)
fig, ax, hm = viz(data_mesh;showfacets=true, color=abs.(data_viz/Z₀))
wgl.Colorbar(fig[1,2],label="|p/Z₀|"); fig

This page was generated using Literate.jl.